Diversity Analysis

Import libraries and data

library(phyloseq)
library(ape)
library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.3.2
## Warning: package 'tidyr' was built under R version 4.3.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::where()  masks ape::where()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(picante)
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 4.3.2
## This is vegan 2.6-4
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
library(ARTool)
library(ggsignif)
library(pairwiseAdonis)
## Loading required package: cluster
library(ggpubr)
## 
## Attaching package: 'ggpubr'
## 
## The following object is masked from 'package:ape':
## 
##     rotate
starch <- read_rds("../21_phyloseq/starch_rarefied_phyloseq.RDS")

Separate human and mouse for the independent analysis

human <- subset_samples(starch, cohort == "human")
mouse <- subset_samples(starch, cohort == "mouse")
pal <- c("#DD4124FF", "#faa478", "#6ac1de", "#0269a1FF")
pal2 <- c("#00346EFF", "#007CBFFF", "#06ABDFFF", "#EDD03EFF", "#F5A200FF", "#6D8600FF", "#424D0CFF")
pal3 <- c("#FDAE61FF", "#ABDDA4FF")

Alpha diversity

richness <- estimate_richness(starch)
rownames(richness) <- rownames(starch@sam_data)
starch_rich <- merge(starch@sam_data, richness, by = 0) %>% dplyr::rename(sample_id = Row.names)

WIP:

otu <- otu_table(starch) %>% as.matrix() 
data_evenness <- diversity(otu) / log(specnumber(otu))

Shannon and Chao1

Statistical Analysis

Merged
shannon_kw <- kruskal.test(Shannon ~ group, data = starch_rich)
shannon_kw
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Shannon by group
## Kruskal-Wallis chi-squared = 26.916, df = 3, p-value = 6.132e-06
if (shannon_kw$p.value < 0.05){
    pairwise.wilcox.test(starch_rich$Shannon, starch_rich$group, 
                       p.adjust = "bonferroni")
}
## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  starch_rich$Shannon and starch_rich$group 
## 
##          human-C human-RS mouse-C
## human-RS 0.058   -        -      
## mouse-C  6.1e-06 0.019    -      
## mouse-RS 9.8e-05 0.218    1.000  
## 
## P value adjustment method: bonferroni

There is a significant difference (p-val < 0.05) in the Shannon index between: - mouse-C & human-C - mouse-C & human-RS - mouse-RS & human-C

shannon_kw <- kruskal.test(Shannon ~ gender, data = starch_rich)
shannon_kw
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Shannon by gender
## Kruskal-Wallis chi-squared = 5.6784, df = 1, p-value = 0.01718

Overall shows a significant difference in the Shannon index between genders.

chao1_kw <- kruskal.test(Chao1 ~ group, data = starch_rich)
chao1_kw
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Chao1 by group
## Kruskal-Wallis chi-squared = 10.757, df = 3, p-value = 0.01311
if (chao1_kw$p.value < 0.05){
    pairwise.wilcox.test(starch_rich$Chao1, starch_rich$group, 
                       p.adjust = "bonferroni")
}
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  starch_rich$Chao1 and starch_rich$group 
## 
##          human-C human-RS mouse-C
## human-RS 1.000   -        -      
## mouse-C  0.018   0.770    -      
## mouse-RS 0.045   1.000    1.000  
## 
## P value adjustment method: bonferroni

Difference in Chao1 index between: human-c & mouse-c and human-C & mouse-rs No difference between: human-RS and mouse-RS!!

chao1_kw <- kruskal.test(Chao1 ~ gender, data = starch_rich)
chao1_kw
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Chao1 by gender
## Kruskal-Wallis chi-squared = 3.33, df = 1, p-value = 0.06803
if (chao1_kw$p.value < 0.05){
    pairwise.wilcox.test(starch_rich$Chao1, starch_rich$gender, 
                       p.adjust = "bonferroni")
}

No significant difference between genders.

Mouse
starch_rich %>% 
    filter(cohort == "mouse") %>%
    kruskal.test(Shannon ~ trt, data = .)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Shannon by trt
## Kruskal-Wallis chi-squared = 1.2027, df = 1, p-value = 0.2728

No significant difference in Shannon index between treatments.

starch_rich %>% 
    filter(cohort == "mouse") %>%
    kruskal.test(Chao1 ~ trt, data = .)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Chao1 by trt
## Kruskal-Wallis chi-squared = 0.98674, df = 1, p-value = 0.3205

No significant difference in Chao1 index between treatments.

starch_rich %>% 
    filter(cohort == "mouse") %>%
    kruskal.test(Shannon ~ gender, data = .)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Shannon by gender
## Kruskal-Wallis chi-squared = 8.3409, df = 1, p-value = 0.003876

Significant difference in Shannon index between genders.

starch_rich %>% 
    filter(cohort == "mouse") %>%
    kruskal.test(Chao1 ~ gender, data = .)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Chao1 by gender
## Kruskal-Wallis chi-squared = 5.6614, df = 1, p-value = 0.01734

Significant difference in Chao1 index between genders.

starch_rich %>% 
    filter(cohort == "mouse") %>%
    kruskal.test(Shannon ~ starch_type, data = .)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Shannon by starch_type
## Kruskal-Wallis chi-squared = 8.2792, df = 4, p-value = 0.08187

No significant difference in Shannon index between starch types

mouse_kw <- starch_rich %>% 
    filter(cohort == "mouse") %>%
    kruskal.test(Chao1 ~ starch_type, data = .)
mouse_rich <- starch_rich %>% 
    filter(cohort == "mouse")

chao1_kw <- mouse_rich %>%
    kruskal.test(Chao1 ~ starch_type, data = .)

if (chao1_kw$p.value < 0.05){
    pairwise.wilcox.test(mouse_rich$Chao1, mouse_rich$starch_type, 
                       p.adjust = "bonferroni")
}
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  mouse_rich$Chao1 and mouse_rich$starch_type 
## 
##     BEP   CKP   CTL   LEN  
## CKP 0.068 -     -     -    
## CTL 1.000 0.092 -     -    
## LEN 1.000 0.459 1.000 -    
## PTB 1.000 0.028 1.000 1.000
## 
## P value adjustment method: bonferroni

Significant difference in mean Chao1 index between CKP and PTB.

Human
starch_rich %>% 
    filter(cohort == "human") %>%
    kruskal.test(Shannon ~ trt, data = .)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Shannon by trt
## Kruskal-Wallis chi-squared = 6.5108, df = 1, p-value = 0.01072

Significant difference in Shannon index between treatments

starch_rich %>% 
    filter(cohort == "human") %>%
    kruskal.test(Chao1 ~ trt, data = .)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Chao1 by trt
## Kruskal-Wallis chi-squared = 1.5583, df = 1, p-value = 0.2119

No significant difference in Chao1 index between treatments.

starch_rich %>% 
    filter(cohort == "human") %>%
    kruskal.test(Shannon ~ gender, data = .)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Shannon by gender
## Kruskal-Wallis chi-squared = 4.6282, df = 1, p-value = 0.03145

Significant difference in Shannon index between genders.

starch_rich %>% 
    filter(cohort == "human") %>%
    kruskal.test(Chao1 ~ gender, data = .)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Chao1 by gender
## Kruskal-Wallis chi-squared = 1.3475, df = 1, p-value = 0.2457

No significant difference in Chao1 index between genders.

Plots

my_comparisons <- list(c("human-C", "human-RS"),
                                             c("mouse-C", "mouse-RS"),
                                             c("mouse-RS","human-RS"),
                                             c("mouse-C","human-C"))
starch_comparisons <- list(c("CTL", "CKP"),
                                                     c("CKP", "LEN"),
                                                     c("CKP", "BEP"),
                                                     c("CKP", "PTB")
                                                     )
Merged
w <- 2400
h <- 2000
starch_rich %>%
    ggboxplot(x = "group", y = "Shannon", color =  "group",
                        palette = pal) +
    stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none",
        axis.title.x = element_blank()
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_group.png", units = "px", width = w, height = h, create.dir = TRUE)

starch_rich %>%
    ggboxplot(x = "gender", y = "Shannon", color =  "gender",
                        palette = pal3) +
    stat_compare_means(label = "p.signif", label.y = 4.2, label.x = 1.5) +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none",
        axis.title.x = element_blank()
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_gender.png", units = "px", width = w, height = h, create.dir = TRUE)

starch_rich %>%
    ggboxplot(x = "gender", y = "Shannon", color =  "gender",
                        palette = pal3, facet.by = "group") +
    stat_compare_means(label = "p.signif", label.y = 4.2, label.x = 1.5) +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none",
        axis.title.x = element_blank()
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_gender_facet.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
    ggboxplot(x = "group", y = "Chao1", color =  "group",
                        palette = pal) +
    stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none",
        axis.title.x = element_blank()
    )
## Warning in wilcox.test.default(c(94.3333333333333, 130.5, 83, 124, 85, 97, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(67, 75, 58, 71, 69.4285714285714, 72, 69.5, :
## cannot compute exact p-value with ties

ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_group.png", units = "px", width = w, height = h, create.dir = TRUE)
## Warning in wilcox.test.default(c(94.3333333333333, 130.5, 83, 124, 85, 97, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(94.3333333333333, 130.5, 83, 124, 85, 97, :
## cannot compute exact p-value with ties
starch_rich %>%
    ggboxplot(x = "gender", y = "Chao1", color =  "gender",
                        palette = pal3) +
    stat_compare_means(label = "p.signif", label.y = 4.2, label.x = 1.5) +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none",
        axis.title.x = element_blank()
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_gender.png", units = "px", width = w, height = h, create.dir = TRUE)

starch_rich %>%
    ggboxplot(x = "gender", y = "Chao1", color =  "gender",
                        palette = pal3, facet.by = "group") +
    stat_compare_means(label = "p.signif", label.y = 4.2, label.x = 1.5) +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none",
        axis.title.x = element_blank()
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_gender_facet.png", units = "px", width = w, height = h, create.dir = TRUE)
Mouse
starch_rich %>%
    filter(cohort == "mouse") %>%
    ggboxplot(x = "gender", y = "Shannon", color =  "gender",
                        palette = pal3) +
    facet_wrap(.~trt) +
    stat_compare_means(method = "kruskal.test",
                                         hide.ns = TRUE,
                                         label = "p.signif") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "right"
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_gender_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

starch_rich %>%
    filter(cohort == "mouse") %>%
    ggboxplot(x = "trt", y = "Shannon", color =  "trt",
                        palette = rev(pal[3:4])) +
    facet_wrap(.~gender) +
    stat_compare_means(method = "kruskal.test",
                                         hide.ns = TRUE,
                                         label = "p.signif") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "right"
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_trt_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

starch_rich %>%
    filter(cohort == "mouse") %>%
    ggboxplot(x = "starch_type", y = "Shannon", color =  "starch_type",
                        palette = pal2) +
    stat_compare_means(method = "kruskal.test",
                                         hide.ns = TRUE,
                                         label = "p.signif") +
    xlab("Starch Type") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none"
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_type_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
    filter(cohort == "mouse") %>%
    ggboxplot(x = "gender", y = "Chao1", color =  "gender",
                        palette = pal3) +
    facet_wrap(.~trt) +
    stat_compare_means(method = "kruskal.test",
                                         hide.ns = TRUE,
                                         label = "p.signif") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "right"
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_gender_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

starch_rich %>%
    filter(cohort == "mouse") %>%
    ggboxplot(x = "trt", y = "Chao1", color =  "trt",
                        palette = rev(pal[3:4])) +
    facet_wrap(.~gender) +
    stat_compare_means(method = "kruskal.test",
                                         hide.ns = TRUE,
                                         label = "p.signif") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "right"
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_trt_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

starch_rich %>%
    filter(cohort == "mouse") %>%
    ggboxplot(x = "starch_type", y = "Chao1", color =  "starch_type",
                        palette = pal2) +
    stat_compare_means(comparisons = starch_comparisons,
                                         method = "kruskal.test",
                                         hide.ns = TRUE,
                                         label = "p.signif") +
    xlab("Starch Type") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none"
    )
## Warning: Computation failed in `stat_signif()`.
## Caused by error in `kruskal.test.default()`:
## ! 'x' and 'g' must have the same length

ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_type_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
## Warning: Computation failed in `stat_signif()`.
## Caused by error in `kruskal.test.default()`:
## ! 'x' and 'g' must have the same length
Human
starch_rich %>%
    filter(cohort == "human") %>%
    ggboxplot(x = "gender", y = "Shannon", color =  "gender",
                        palette = pal3) +
    facet_wrap(.~trt) +
    stat_compare_means(method = "kruskal.test",
                                         hide.ns = TRUE,
                                         label = "p.signif") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "right"
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_gender_human.png", units = "px", width = w, height = h, create.dir = TRUE)

starch_rich %>%
    filter(cohort == "human") %>%
    ggboxplot(x = "trt", y = "Shannon", color =  "trt",
                        palette = pal) +
    facet_wrap(.~gender) +
    stat_compare_means(method = "kruskal.test", 
                                         hide.ns = TRUE,
                                         label = "p.signif") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "right"
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_trt_human.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
    filter(cohort == "human") %>%
    ggboxplot(x = "gender", y = "Chao1", color =  "gender",
                        palette = pal3) +
    facet_wrap(.~trt) +
    stat_compare_means(method = "kruskal.test",
                                         hide.ns = TRUE,
                                         label = "p.signif") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "right"
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_gender_human.png", units = "px", width = w, height = h, create.dir = TRUE)

starch_rich %>%
    filter(cohort == "human") %>%
    ggboxplot(x = "trt", y = "Chao1", color =  "trt",
                        palette = pal) +
    facet_wrap(.~gender) +
    stat_compare_means(method = "kruskal.test", 
                                         hide.ns = TRUE,
                                         label = "p.signif") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "right"
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_trt_human.png", units = "px", width = w, height = h, create.dir = TRUE)

Faith’s Phylogenetic Diversity

phylo_dist <- pd(t(otu_table(starch)), phy_tree(starch),
                 include.root=F) 

# add PD to metadata table
starch_rich$PD <- phylo_dist$PD

Statistical Analysis

Group
pd_kw <- kruskal.test(PD ~ group, data = starch_rich)
pd_kw
## 
##  Kruskal-Wallis rank sum test
## 
## data:  PD by group
## Kruskal-Wallis chi-squared = 4.5904, df = 3, p-value = 0.2044
if (pd_kw$p.value < 0.05){
    pairwise.wilcox.test(starch_rich$PD, starch_rich$group, 
                       p.adjust = "bonferroni")
}

There is no significant difference in Faith’s phylogenetic distance between the groups.

Gender
kruskal.test(PD ~ gender, data = starch_rich)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  PD by gender
## Kruskal-Wallis chi-squared = 1.2316, df = 1, p-value = 0.2671

No difference between genders.

Mouse
starch_rich %>% 
    filter(cohort == "mouse") %>%
    kruskal.test(PD ~ gender, data = .)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  PD by gender
## Kruskal-Wallis chi-squared = 3.1929, df = 1, p-value = 0.07396

No difference between genders.

starch_rich %>% 
    filter(cohort == "mouse") %>%
    kruskal.test(PD ~ starch_type, data = .)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  PD by starch_type
## Kruskal-Wallis chi-squared = 6.8053, df = 4, p-value = 0.1465

No difference between starch types.

Human
starch_rich %>% 
    filter(cohort == "human") %>%
    kruskal.test(PD ~ gender, data = .)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  PD by gender
## Kruskal-Wallis chi-squared = 1.0385, df = 1, p-value = 0.3082

No difference between genders.

Plots

Merged
starch_rich %>%
    ggboxplot(x = "group", y = "PD", color =  "group",
                        palette = pal) +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none",
        axis.title.x = element_blank()
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/pd_group.png", units = "px", width = w, height = h, create.dir = TRUE)

starch_rich %>%
    ggboxplot(x = "gender", y = "PD", color =  "gender",
                        palette = pal3) +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none"
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/pd_gender.png", units = "px", width = w, height = h, create.dir = TRUE)

starch_rich %>%
    ggboxplot(x = "gender", y = "PD", color =  "gender",
                        palette = pal3, facet.by = "group") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none"
    )

ggsave("../30_diversity-analysis/plots/alpha-diversity/pd_gender_facet.png", units = "px", width = w, height = h, create.dir = TRUE)

Beta Diversity

starchdf <- as.data.frame(sample_data(starch))
humandf <- as.data.frame(sample_data(human))
mousedf <- as.data.frame(sample_data(mouse))

Bray Curtis

Merged

bc_dm <- phyloseq::distance(starch, method="bray")
adonis2(bc_dm ~ group*gender, data = data.frame(starchdf))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = bc_dm ~ group * gender, data = data.frame(starchdf))
##               Df SumOfSqs      R2      F Pr(>F)    
## group          3    9.325 0.25341 13.304  0.001 ***
## gender         1    3.421 0.09296 14.641  0.001 ***
## group:gender   3    2.323 0.06313  3.314  0.001 ***
## Residual      93   21.729 0.59050                  
## Total        100   36.797 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(bc_dm ~ trt*cohort*gender, data = data.frame(starchdf))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = bc_dm ~ trt * cohort * gender, data = data.frame(starchdf))
##                    Df SumOfSqs      R2       F Pr(>F)    
## trt                 1    1.737 0.04722  7.4361  0.001 ***
## cohort              1    7.253 0.19710 31.0412  0.001 ***
## gender              1    3.451 0.09379 14.7715  0.001 ***
## trt:cohort          1    0.304 0.00827  1.3028  0.208    
## trt:gender          1    0.983 0.02670  4.2055  0.001 ***
## cohort:gender       1    0.956 0.02597  4.0898  0.001 ***
## trt:cohort:gender   1    0.385 0.01046  1.6468  0.109    
## Residual           93   21.729 0.59050                   
## Total             100   36.797 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There is a significant difference between treatments, cohorts, genders, at least two trt:gender pair, and at least two cohort:gender pair.

pcoa_bc <- ordinate(starch, method="PCoA", distance=bc_dm)

plot_ordination(starch, pcoa_bc, color = "group") +
    scale_color_manual(values = pal) +
  labs(col = "Group") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/bc_group.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(starch, pcoa_bc, color = "gender") +
    scale_color_manual(values = pal3) +
  labs(col = "Gender") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/bc_gender.png", units = "px", width = w, height = h, create.dir = TRUE)

Human

human_bc_dm <- distance(human, method="bray")
adonis2(human_bc_dm ~ trt*gender, data = data.frame(humandf))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = human_bc_dm ~ trt * gender, data = data.frame(humandf))
##            Df SumOfSqs      R2      F Pr(>F)  
## trt         1   0.2863 0.04489 1.0845  0.322  
## gender      1   0.4024 0.06310 1.5246  0.028 *
## trt:gender  1   0.1455 0.02281 0.5512  0.993  
## Residual   21   5.5430 0.86919                
## Total      24   6.3772 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There is only a significant difference in community structure between genders. I should probably consider here the effect of trtseq (fixed effect) or individuals as random error….

human_pcoa_bc <- ordinate(human, method="PCoA", distance=human_bc_dm)

plot_ordination(human, human_pcoa_bc, color = "trt") +
    scale_color_manual(values = pal) +
  labs(col = "Treatment") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/bc_trt_human.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(human, human_pcoa_bc, color = "gender") +
    scale_color_manual(values = pal3) +
  labs(col = "Gender") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/bc_gender_human.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(human, human_pcoa_bc, color = "gender") +
    scale_color_manual(values = pal3) +
    facet_grid(.~trt) +
  labs(col = "Gender") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/bc_gender_facet_human.png", units = "px", width = w, height = h, create.dir = TRUE)

Mouse

mouse_bc_dm <- distance(mouse, method="bray")
adonis2(mouse_bc_dm ~ trt*starch_type*gender, data = data.frame(mousedf))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = mouse_bc_dm ~ trt * starch_type * gender, data = data.frame(mousedf))
##                    Df SumOfSqs      R2       F Pr(>F)    
## trt                 1   0.4224 0.01937  2.6609  0.030 *  
## starch_type         3   3.0381 0.13934  6.3789  0.001 ***
## gender              1   4.2557 0.19518 26.8060  0.001 ***
## trt:gender          1   0.8667 0.03975  5.4595  0.001 ***
## starch_type:gender  3   2.7427 0.12579  5.7587  0.001 ***
## Residual           66  10.4781 0.48056                   
## Total              75  21.8038 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pw <- pairwiseAdonis::pairwise.adonis2(mouse_bc_dm ~ trt*starch_type*gender, 
                                                                             data = data.frame(mousedf))

print("significantly different pairs: ")
## [1] "significantly different pairs: "
for (i in range(2,(length(pw)))){
    pair <- pw[[i]]
    if (pair$`Pr(>F)`[1] < 0.05) {
        print(names(pw)[i])
    }
}
## [1] "RS_vs_C"
## [1] "RS_vs_C"

There is a significant difference in community structure between treatments, at least 2 starch types, genders, at least two trt:gender pairs, and at least two starch_type:gender pairs.

mouse_pcoa_bc <- ordinate(mouse, method="PCoA", distance=mouse_bc_dm)

plot_ordination(mouse, mouse_pcoa_bc, color = "trt") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal[3:4]) +
  labs(col = "Treatment") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/bc_trt_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(mouse, mouse_pcoa_bc, color = "starch_type") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal2) +
  labs(col = "Starch Type") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/bc_type_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(mouse, mouse_pcoa_bc, color = "gender") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal3) +
  labs(col = "Starch Type") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/bc_gender_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

Jaccard

Merged

jc_dm <- phyloseq::distance(starch, method = "jaccard", binary = TRUE)
adonis2(jc_dm ~ group*gender, data = data.frame(starchdf))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = jc_dm ~ group * gender, data = data.frame(starchdf))
##               Df SumOfSqs      R2      F Pr(>F)    
## group          3    7.162 0.17428 7.0604  0.001 ***
## gender         1    1.194 0.02906 3.5315  0.001 ***
## group:gender   3    1.292 0.03144 1.2736  0.042 *  
## Residual      93   31.444 0.76522                  
## Total        100   41.092 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(jc_dm ~ cohort*trt*gender, data = data.frame(starchdf))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = jc_dm ~ cohort * trt * gender, data = data.frame(starchdf))
##                    Df SumOfSqs      R2       F Pr(>F)    
## cohort              1    6.352 0.15458 18.7866  0.001 ***
## trt                 1    0.439 0.01070  1.2999  0.094 .  
## gender              1    1.204 0.02930  3.5609  0.001 ***
## cohort:trt          1    0.360 0.00877  1.0655  0.260    
## cohort:gender       1    0.641 0.01559  1.8947  0.010 ** 
## trt:gender          1    0.375 0.00912  1.1082  0.237    
## cohort:trt:gender   1    0.277 0.00673  0.8180  0.839    
## Residual           93   31.444 0.76522                   
## Total             100   41.092 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There is a significant difference in community composition between cohorts, genders, and at least 2 cohort:gender pairs.

pcoa_jc <- ordinate(starch, method="PCoA", distance=jc_dm)

plot_ordination(starch, pcoa_jc, color = "group") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal) +
  labs(col = "Group") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/jc_group.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(starch, pcoa_jc, color = "gender") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal3) +
  labs(col = "Gender") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/jc_gender.png", units = "px", width = w, height = h, create.dir = TRUE)

Human

human_jc_dm <- distance(human, method = "jaccard", binary = TRUE)
adonis2(human_jc_dm ~ trt*gender, data = data.frame(humandf))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = human_jc_dm ~ trt * gender, data = data.frame(humandf))
##            Df SumOfSqs      R2      F Pr(>F)  
## trt         1   0.2831 0.03702 0.8769  0.762  
## gender      1   0.4030 0.05271 1.2484  0.092 .
## trt:gender  1   0.1805 0.02361 0.5592  1.000  
## Residual   21   6.7794 0.88666                
## Total      24   7.6461 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Community composition is not significantly different between treatments or gender.

human_pcoa_jc <- ordinate(human, method="PCoA", distance=human_jc_dm)

plot_ordination(human, human_pcoa_jc, color = "trt") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal) +
  labs(col = "Treatment") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/jc_trt_human.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(human, human_pcoa_jc, color = "gender") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal3) +
  labs(col = "Gender") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/jc_gender_human.png", units = "px", width = w, height = h, create.dir = TRUE)

Mouse

mouse_jc_dm <- distance(mouse, method = "jaccard", binary = TRUE)
adonis2(mouse_jc_dm ~ gender*trt*starch_type, data = data.frame(mousedf))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = mouse_jc_dm ~ gender * trt * starch_type, data = data.frame(mousedf))
##                    Df SumOfSqs      R2      F Pr(>F)    
## gender              1   1.4456 0.05335 4.4497  0.001 ***
## trt                 1   0.5127 0.01892 1.5781  0.004 ** 
## starch_type         3   1.7031 0.06286 1.7475  0.001 ***
## gender:trt          1   0.4717 0.01741 1.4518  0.005 ** 
## gender:starch_type  3   1.5197 0.05609 1.5593  0.001 ***
## Residual           66  21.4413 0.79137                  
## Total              75  27.0939 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pw <- pairwiseAdonis::pairwise.adonis2(mouse_jc_dm ~ starch_type, data = data.frame(mousedf))

print("significantly different pairs: ")
## [1] "significantly different pairs: "
for (i in range(2,(length(pw)))){
    pair <- pw[[i]]
    if (pair$`Pr(>F)`[1] < 0.05) {
        print(names(pw)[i])
    }
}
## [1] "PTB_vs_BEP"
## [1] "CKP_vs_CTL"

Community composition is significantly different between genders, treatments, starch types, at least 2 gender:trt pairs and at least 2 gender:starch_type pairs. The significantly different starch_type groups are PTB-BEP and CKP-CTL.

mouse_pcoa_jc <- ordinate(mouse, method="PCoA", distance=mouse_jc_dm)

plot_ordination(mouse, mouse_pcoa_jc, color = "trt") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal[3:4]) +
  labs(col = "Treatment") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/jc_trt_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(mouse, mouse_pcoa_jc, color = "starch_type") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal2) +
  labs(col = "Starch Type") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/jc_type_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(mouse, mouse_pcoa_jc, color = "gender") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal3) +
  labs(col = "Gender") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/jc_gender_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

Unweighted Unifrac

Merged

uu_dm <- phyloseq::distance(starch, method = "unifrac")
## Warning in matrix(tree$edge[order(tree$edge[, 1]), ][, 2], byrow = TRUE, : data
## length [3823] is not a sub-multiple or multiple of the number of rows [1912]
adonis2(uu_dm ~ group*gender, data = data.frame(starchdf))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = uu_dm ~ group * gender, data = data.frame(starchdf))
##               Df SumOfSqs      R2      F Pr(>F)    
## group          3   4.1175 0.21931 9.5810  0.001 ***
## gender         1   0.6622 0.03527 4.6226  0.001 ***
## group:gender   3   0.6730 0.03585 1.5660  0.031 *  
## Residual      93  13.3224 0.70958                  
## Total        100  18.7751 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(uu_dm ~ cohort*trt*gender, data = data.frame(starchdf))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = uu_dm ~ cohort * trt * gender, data = data.frame(starchdf))
##                    Df SumOfSqs      R2       F Pr(>F)    
## cohort              1   3.8031 0.20256 26.5482  0.001 ***
## trt                 1   0.1644 0.00876  1.1478  0.282    
## gender              1   0.6630 0.03532  4.6286  0.001 ***
## cohort:trt          1   0.1491 0.00794  1.0409  0.356    
## cohort:gender       1   0.3743 0.01994  2.6129  0.008 ** 
## trt:gender          1   0.1801 0.00959  1.2573  0.225    
## cohort:trt:gender   1   0.1186 0.00632  0.8278  0.642    
## Residual           93  13.3224 0.70958                   
## Total             100  18.7751 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There is a significant difference in community composition and phylogeny between cohorts, genders and at least 2 group:gender pairs.

pcoa_uu <- ordinate(starch, method="PCoA", distance=uu_dm)

plot_ordination(starch, pcoa_uu, color = "group") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal) +
  labs(col = "Group") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/uu_group.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(starch, pcoa_uu, color = "gender") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal3) +
  labs(col = "Gender") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/uu_gender.png", units = "px", width = w, height = h, create.dir = TRUE)

Human

human_uu_dm <- distance(human, method = "unifrac")
## Warning in matrix(tree$edge[order(tree$edge[, 1]), ][, 2], byrow = TRUE, : data
## length [3823] is not a sub-multiple or multiple of the number of rows [1912]
adonis2(human_uu_dm ~ trt*gender, data = data.frame(humandf), permutations = 9999)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = human_uu_dm ~ trt * gender, data = data.frame(humandf), permutations = 9999)
##            Df SumOfSqs      R2      F Pr(>F)  
## trt         1   0.1656 0.04265 1.0377 0.3880  
## gender      1   0.2459 0.06336 1.5414 0.0456 *
## trt:gender  1   0.1196 0.03081 0.7497 0.8319  
## Residual   21   3.3503 0.86318                
## Total      24   3.8814 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There is a significant difference in community composition and phylogeny between genders.

human_pcoa_uu <- ordinate(human, method="PCoA", distance=human_uu_dm)

plot_ordination(human, human_pcoa_uu, color = "trt") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal) +
  labs(col = "Treatment") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/uu_trt_human.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(human, human_pcoa_uu, color = "gender") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal3) +
  labs(col = "Gender") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/uu_gender_human.png", units = "px", width = w, height = h, create.dir = TRUE)

Mouse

mouse_uu_dm <- distance(mouse, method = "unifrac")
## Warning in matrix(tree$edge[order(tree$edge[, 1]), ][, 2], byrow = TRUE, : data
## length [3823] is not a sub-multiple or multiple of the number of rows [1912]
adonis2(mouse_uu_dm ~ trt*gender*starch_type, data = data.frame(mousedf))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = mouse_uu_dm ~ trt * gender * starch_type, data = data.frame(mousedf))
##                    Df SumOfSqs      R2      F Pr(>F)    
## trt                 1   0.1488 0.01342 1.1175  0.301    
## gender              1   0.7906 0.07128 5.9357  0.001 ***
## starch_type         3   0.5704 0.05143 1.4276  0.075 .  
## trt:gender          1   0.1799 0.01622 1.3506  0.164    
## gender:starch_type  3   0.6101 0.05501 1.5269  0.030 *  
## Residual           66   8.7907 0.79263                  
## Total              75  11.0906 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There no significant difference in community composition and phylogeny between treatments but there is a difference between genders and at least 2 gender:starch_type pairs.

mouse_pcoa_uu <- ordinate(mouse, method="PCoA", distance=mouse_uu_dm)

plot_ordination(mouse, mouse_pcoa_uu, color = "trt") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal[3:4]) +
  labs(col = "Treatment") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/uu_trt_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(mouse, mouse_pcoa_uu, color = "starch_type") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal2) +
  labs(col = "Starch Type") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/uu_type_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(mouse, mouse_pcoa_uu, color = "gender") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal3) +
  labs(col = "Gender") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/uu_gender_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

there’s an interesting clustering, but we’d need more metadata to see what is causing this difference.

Weighted Unifrac

Merged

wu_dm <- phyloseq::distance(starch, method = "wunifrac")
## Warning in matrix(tree$edge[order(tree$edge[, 1]), ][, 2], byrow = TRUE, : data
## length [3823] is not a sub-multiple or multiple of the number of rows [1912]
adonis2(wu_dm ~ group*gender, data = data.frame(starchdf), permutations = 9999)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = wu_dm ~ group * gender, data = data.frame(starchdf), permutations = 9999)
##               Df SumOfSqs      R2      F Pr(>F)    
## group          3   4.0276 0.44175 30.339 0.0001 ***
## gender         1   0.5911 0.06483 13.357 0.0001 ***
## group:gender   3   0.3834 0.04205  2.888 0.0084 ** 
## Residual      93   4.1153 0.45137                  
## Total        100   9.1173 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(wu_dm ~ cohort*trt*gender, data = data.frame(starchdf), permutations = 9999)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = wu_dm ~ cohort * trt * gender, data = data.frame(starchdf), permutations = 9999)
##                    Df SumOfSqs      R2       F Pr(>F)    
## cohort              1   3.9009 0.42786 88.1552 0.0001 ***
## trt                 1   0.0558 0.00612  1.2619 0.2566    
## gender              1   0.5791 0.06352 13.0872 0.0001 ***
## cohort:trt          1   0.0828 0.00908  1.8705 0.1291    
## cohort:gender       1   0.1137 0.01247  2.5698 0.0636 .  
## trt:gender          1   0.2141 0.02348  4.8382 0.0073 ** 
## cohort:trt:gender   1   0.0556 0.00610  1.2559 0.2478    
## Residual           93   4.1153 0.45137                   
## Total             100   9.1173 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Significant difference in community structure and phylogeny between cohorts, genders and at least 2 trt:gender pairs.

pcoa_wu <- ordinate(starch, method="PCoA", distance=wu_dm)

plot_ordination(starch, pcoa_wu, color = "group") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal) +
  labs(col = "Group") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/wu_group.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(starch, pcoa_wu, color = "gender") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal3) +
  labs(col = "Gender") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/wu_gender.png", units = "px", width = w, height = h, create.dir = TRUE)

Human

human_wu_dm <- distance(human, method = "wunifrac")
## Warning in matrix(tree$edge[order(tree$edge[, 1]), ][, 2], byrow = TRUE, : data
## length [3823] is not a sub-multiple or multiple of the number of rows [1912]
adonis2(human_wu_dm ~ trt*gender, data = data.frame(humandf), permutations = 9999)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = human_wu_dm ~ trt * gender, data = data.frame(humandf), permutations = 9999)
##            Df SumOfSqs      R2      F Pr(>F)
## trt         1  0.02845 0.04848 1.2105 0.2718
## gender      1  0.03558 0.06063 1.5141 0.1968
## trt:gender  1  0.02929 0.04992 1.2465 0.2612
## Residual   21  0.49351 0.84097              
## Total      24  0.58683 1.00000

There is no significant difference in community structure and phylogeny between treatments, genders or the interaction between them.

human_pcoa_wu <- ordinate(human, method="PCoA", distance=human_wu_dm)

plot_ordination(human, human_pcoa_wu, color = "trt") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal[1:2]) +
  labs(col = "Treatment") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/wu_trt_human.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(human, human_pcoa_wu, color = "gender") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal3) +
  labs(col = "Gender") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/wu_gender_human.png", units = "px", width = w, height = h, create.dir = TRUE)

Mouse

mouse_wu_dm <- distance(mouse, method = "wunifrac")
## Warning in matrix(tree$edge[order(tree$edge[, 1]), ][, 2], byrow = TRUE, : data
## length [3823] is not a sub-multiple or multiple of the number of rows [1912]
adonis2(mouse_wu_dm ~ trt*gender*starch_type, data = data.frame(mousedf))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = mouse_wu_dm ~ trt * gender * starch_type, data = data.frame(mousedf))
##                    Df SumOfSqs      R2       F Pr(>F)    
## trt                 1   0.0982 0.02121  2.7985  0.027 *  
## gender              1   0.6692 0.14455 19.0702  0.001 ***
## starch_type         3   0.6250 0.13501  5.9370  0.001 ***
## trt:gender          1   0.2444 0.05279  6.9650  0.001 ***
## gender:starch_type  3   0.6767 0.14616  6.4276  0.001 ***
## Residual           66   2.3161 0.50028                   
## Total              75   4.6296 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pw <- pairwiseAdonis::pairwise.adonis2(mouse_wu_dm ~ starch_type, data = data.frame(mousedf))

print("significantly different pairs: ")
## [1] "significantly different pairs: "
for (i in range(2,(length(pw)))){
    pair <- pw[[i]]
    if (pair$`Pr(>F)`[1] < 0.05) {
        print(names(pw)[i])
    }
}
## [1] "PTB_vs_BEP"
## [1] "CKP_vs_CTL"

There is a significant difference in community structure and phylogeny between treatments, genders, at least 2 starch types, at least 2 trt:gender pairs, and at least 2 gender:starch_type pairs. The significantly different starch type pairs are: PTB-BEP and CKP-CTL.

mouse_pcoa_wu <- ordinate(mouse, method="PCoA", distance=mouse_wu_dm)

plot_ordination(mouse, mouse_pcoa_wu, color = "trt") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal[3:4]) +
  labs(col = "Treatment") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/wu_trt_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(mouse, mouse_pcoa_wu, color = "starch_type") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal2) +
  labs(col = "Starch Type") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/wu_type_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)

plot_ordination(mouse, mouse_pcoa_wu, color = "gender") +
    #stat_ellipse(alpha = 0.5) +
    scale_color_manual(values = pal3) +
  labs(col = "Gender") +
    theme_bw() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    )

ggsave("../30_diversity-analysis/plots/beta-diversity/wu_gender_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)